function explicit_sin2

%  Solves the heat equation using explicit method
%       diff(u,x,x) = diff(u,t) + f(x,t)   for xL < x < xr, 0 < t < tmax
%  where
%      u = 0  at x=xL,xR  and  u = g(x) at t = 0

%  g(x)=sin(2*pi*x)  and f(x,t)=0

% clear all previous variables and plots
clear *
clf

% set parameters
N=20;
M=90;
tmax=0.1;
xL=0;
xR=1;

fprintf('\n    Solution Computed with N = %3.0f and M = %4.0f\n\n',N,M)

% generate the points along the x-axis, x(1)=xL and x(N+2)=xR
x=linspace(xL,xR,N+2);
h=x(2)-x(1);

% generate the points along the t-axis, t(1)=0 and t(M+1)=tmax
t=linspace(0,tmax,M+1);
k=t(2)-t(1);
lamda=k/h^2;

% pick time points for plotting (by picking the index)
itotal=3;
it(1)=1+round(0.02/k);
it(2)=1+round(0.04/k);
it(3)=M;
tt(1)=t(it(1)); tt(2)=t(it(2)); tt(3)=t(it(3));

fprintf('\n    Lamda = %5.2e\n\n',lamda)

% find numerical solution
ue=explicit2(x,t,N+2,M+1,h,k,lamda);

% get(gcf)
%set(gcf,'Position', [1155 497 575 470]);
plotsize(560,725)

% plot results
xx=linspace(xL,xR,100);
hold on
for itt=1:itotal
	
	subplot(3,1,4-itt)
	
	% plot numerical and exact solutions
	plot(x,ue(:,it(itt)),'r')
	hold on
	u=exp(-4*pi*pi*tt(itt))*sin(2*pi*xx);
	plot(xx,u,'k')

	% define axes used in plot
	xlabel('x-axis','FontSize',14,'FontWeight','bold')
	ylabel('Solution','FontSize',14,'FontWeight','bold')
	
	% have MATLAB use certain plot options (all are optional)
	box on
	set(gca,'FontSize',14);
	if itt==1
		say=['Time = 0.02']; 
		legend([' M = ',num2str(M,'%3.0f')],' Exact',3);
		set(findobj(gcf,'tag','legend'),'FontSize',12,'FontWeight','bold');
	elseif itt==2
		say=['Time = 0.04']; 
	else
		say=['Time = 0.1']; 
	end;
	ylim;
	yt=0.75*ans(2);
	text(0.75,yt,say,'FontSize',14,'FontWeight','bold')
	hold off
	
end;
say=['Heat Equation: exact vs explicit method when u(x,0)=sin(2\pix)'];
title(say,'FontSize',14,'FontWeight','bold')

% explicit method
function UE=explicit2(x,t,N,M,h,k,lamda)
UE=zeros(N,M);
for i=1:N
	UE(i,1)=g(x(i));
end;
for j=2:M
	for i=2:N-1
	  UE(i,j)=lamda*UE(i+1,j-1)+(1-2*lamda)*UE(i,j-1)+lamda*UE(i-1,j-1)-k*f(x(i),t(j-1));
	end;
end;

% subfunction f(x,t)
function q=f(x,t)
q=0;

% subfunction g(x)
function q=g(x)
q=sin(2*pi*x);

% tridiagonal solver
function y = tridiag( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end;
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end;


% subfunction plotsize
function plotsize(width,height)
siz=get(0,'ScreenSize');
bottom=max(siz(4)-height-95,1);
set(gcf,'Position', [2 bottom width height]);